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The Cosmic Microwave Background (CMB), which permeates the 
entire Universe, is the radiation left over from just 380,000 years 
after the Big Bang. On very large scales, the CMB radiation field 
is smooth and isotropic, but the existence of structure in the 
Universe - stars, galaxies, clusters of galaxies, ... - suggests that 
the field should fluctuate on smaller scales. Recent observations, 
from the Cosmic Microwave Background Explorer to the Wilkin- 
son Microwave Anisotropy Project, have strikingly confirmed this 
prediction. 

CMB fluctuations provide clues to the Universe's structure and 
composition shortly after the Big Bang that are critical for testing 
cosmological models. For example, CMB data can be used to 
determine what portion of the Universe is composed of ordinary 
matter versus the mysterious dark matter and dark energy. To 

this end, cosmologists usually summarize the fluctuations by the 
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power spectrum, which gives the variance as a function of angular 
frequency. The spectrum's shape, and in particular the location 
and height of its peaks, relates directly to the parameters in the 
cosmological models. Thus, a critical statistical question is how 
accurately can these peaks be estimated. 

We use recently developed techniques to construct a nonpara- 
metric confidence set for the unknown CMB spectrum. Our esti- 
mated spectrum, based on minimal assumptions, closely matches 
the model-based estimates used by cosmologists, but we can make 
a wide range of additional inferences. We apply these techniques 
to test various models and to extract confidence intervals on cos- 
mological parameters of interest. Our analysis shows that, even 
without parametric assumptions, the first peak is resolved accu- 
rately with current data but that the second and third peaks are 
not. 

Key words and phrases: Confidence sets, nonparametric regres- 
sion, cosmology. 
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1 Introduction 



The "Big Bang" model is misnamed, as one might expect when a term is 
coined as an insult. Cosmologist Fred Hoyle first used the name in a BBC 
radio interview to denigrate the theory, which opposed the then-dominant 
Steady State model. The name Big Bang stuck, as did its evocation of a 
mighty explosion in space. But the image of an explosion is highly misleading. 
What the model actually posits is that the Universe began hot, dense, and 
expanding. 

Within the first second, roughly 13.7 billion years ago, the Universe 
achieved temperatures on the order of 1 trillion Kelvin (K, degrees above 
absolute zero; Schwarz, 2003). The density during that second was high 
enough to stop neutrinos, which interact so weakly with matter that they 
can pass unmolested through a quadrillion kilometers of lead. What ties this 
hot, dense beginning to the Universe we see today is expansion. A useful 
metaphor for the expanding universe is the surface of an inflating balloon. 
As the balloon inflates, space-time itself is stretched; every point moves away 
from every other point. Density falls as the universe expands. If you picture 
a wave oscillating over the surface of the balloon, the wavelength increases. 
Increasing the wavelength of light corresponds to reducing its temperature. 
The Universe thus cools as it expands. 

Within the first three minutes, the Universe's temperature was over one 
billion K. The energy density in space was so high that atoms could not form. 
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Space was filled with a stew of photons, baryons (e.g., protons and neutrons), 
electrons, neutrinos, and other matter. As the temperature cooled below 
1 billion K, light-element nuclei (deuterium, helium, some lithium) formed 
as well, in proportions that fit well with observations. During this period, 
photons (radiation) were the dominant form of energy in the Universe. Any 
fluctuations in density caused by gravity (which affects light and matter) 
were quickly smoothed out and so could not grow. 

When the temperature of the primordial photons had fallen below ap- 
proximately 12,000 K, photons were no longer dominating the interactions 
among all particles. Photons and baryons became coupled in a mathemat- 
ically perfect fluid, while exotic kinds of matter began to clump under the 
influence of gravity. The interaction between this photon-baryon fluid and 
such gravitational overdensities are of critical importance and will be de- 
scribed below. 

When the temperature reached about 3000 K, roughly 380,000 years after 
the Big Bang, electrons and protons could combine to form atoms. This 
decoupled the photon-baryon fluid, and the photons flew free through space. 
This period is named recombination and happened, in cosmic terms, very 
quickly. After another 200 million years, hydrogen formed after recombination 
had clumped enough for the first stars to form, which began the synthesis of 
heavy elements and the formation of galaxies that we see today. 

Most of the photons released at recombination have travelled through 
space for billions of years without interacting with matter. The temperature 
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of these primordial photons has now cooled to about 2.7K, barely above ab- 
solute zero, which puts them in the microwave part of the electromagnetic 
spectrum. This primordial radiation field, which still pervades the Universe, is 
called the Cosmic Microwave Background (CMB). The CMB thus provides 
a snapshot of the moment of recombination, and fluctuations in the tem- 
perature across the sky contain information about the physics of the early 
universe. 

1.1 The Cosmic Microwave Background Radiation 

As we will explain in the remainder of this section, the temperature fluctu- 
ations in the CMB give a snapshot of the physics in the early Universe and 
provide critical tests of cosmological models. In 1992, the Cosmic Microwave 
Background Explorer (COBE) satellite discovered fluctuations in the black- 
body temperature of the CMB (Smoot et al. 1992). These fluctuations are 
small: approximately one thousandth of the mean temperature over the sky. 
Indeed, almost thirty years of experiments since the CMB's discovery could 
not detect any deviation from uniformity. During the ten years following 
COBE, many more refined measurements were taken; notable experiments 
include MAXIMA, DASI, BOOMERANG (Lee et al. 2002, Halverson et al. 
2002, Netterfield et al. 2002). In 2003, the Wilkinson Microwave Anisotropy 
Project (WMAP) considerably refined the picture, increasing spatial resolu- 
tion by a factor of 33 and sensitivity by a factor of 45 over COBE (Bennett et 
al. 2003). In Figure ^ we compare the COBE and WMAP temperature sky 
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maps after removing the mean temperature T = 2.726 Kelvin and adjusting 
for the motion of our galaxy through the Universe. The fluctuations' mag- 
nitudes are just right to explain the large-scale structure in the Universe we 
see today. For example, if they had been much smaller, there would not be 
enough local concentration of mass to seed the formation of galaxies, galaxy 
clusters, et cetera. 

Perhaps the most important summary of the temperature measurements 
used by cosmologists is the power spectrum, which gives the temperature 
variance as a function of spatial frequency. The spectrum's shape, and in 
particular the location and height of its peaks, relates directly to the pa- 
rameters in cosmological models. (See Appendix 2 for a description of these 
parameters.) Thus, a critical statistical question is how accurately can these 
peaks be estimated. Of particular interest are the height and location of the 
first peak and the relative heights of the successive peaks. 

Figure 121 displays an estimated spectrum commonly used by cosmologists 
and highlights the peaks of interest. We will give a more precise definition of 
the spectrum in Section [Ql but here we want to explain how the spectrum's 
shape relates to the physics in the time up to recombination. 

A key to understanding the physics before recombination is, as men- 
tioned earlier, that photons and baryons became coupled into a (perfect) 
fluid. Mathematical techniques for studying fluid dynamics apply well in this 
scenario and have been investigated by many authors (see for instance Hu and 
Sugiyama 1995; Hu 1999, 2001, 2003; Hu and Dodelson 2002). The properties 
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of the fluid are determined by the relative density of photons and baryons in 
the fluid. Photons provided pressure, and the baryons provided inertia. As 
the fluid falls into a gravitational potential well around a clump of higher 
density, the pressure from the photons resists compression and the inertia of 
the baryons increases it. (Large, isolated potential wells were likely rare in 
the early universe; instead, there were random density fluctuations at many 
scales.) The result is an oscillation that produces pressure waves - sound - in 
the photon-baryon fluid. These acoustic oscillations account for much of the 
interesting structure in the spectrum, particularly the size and arrangement 
of peaks. The imprint of those waves remains in the CMB as a pattern of hot 
and cold spots. 

To understand the peaks in the power spectrum, it is helpful to decom- 
pose the acoustic oscillations into their basic components, or modes. The first 
peak of the spectrum represents the fundamental tone of the oscillations, and 
the other peaks in the spectrum represent harmonics of this tone. The funda- 
mental corresponds to the mode for which one compression occurs between 
the Big Bang and recombination. Each successive harmonic corresponds to 
an additional half-cycle, compressions followed by rarefication (decompres- 
sion). Thus, the second peak represents modes that had time to compress 
and then rarefy before the photons were released from the photon-baryon 
fluid. The third peak represents compression-rarefication-compression, and 
so on. 

The height of the first peak is determined by the total energy density. 
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Roughly, with more matter, the gravitational attraction requires more force 
to counteract, deepening the compression and thus increasing the amplitude 
of oscillation. 

Now suppose we increase the density of baryons in the photon-baryon 
fluid. This increases the inertia of the fluid, deepening each compression phase 
without changing the rarefication. The oscillations become asymmetric. What 
this means is the odd-numbered peaks, whose modes end on a compression, 
are enhanced relative to the even-numbered peaks, whose modes end on a 
rarefication. Thus as the baryon fraction increases we should (over some 
range) see a differential effect on the odd and even numbered peaks. 

The third peak in the spectrum provides the clearest support for the exis- 
tence of "dark matter" - a substance of unknown composition that interacts 
at most weakly with baryons (e.g., neutrons, protons) or with photons (that's 
why its dark). To see why, it is illuminating to compare the oscillations in 
two example cases. In the "radiation-dominated era," when photons were the 
dominant form of interaction in the universe, density fluctuations were short- 
lived and unstable. A compressed region of photon-baryon fluid would rarefy 
as described earlier, but as it did so, the overdensity that caused the original 
gravitational well would disappear. Thus, in this case, at most one cycle of 
oscillation would occur between the Big Bang and recombination. We would 
see only one small peak in the temperature power spectrum corresponding to 
the mode (component of oscillation) that reaches maximum compression at 
the time of recombination. When the matter fraction is low, the peak would 
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be small, increasing with the baryon fraction (inertia). 

In the "matter-dominated era," however, most of the energy density was 
in the form of dark matter. The rarefication phase of the oscillation would 
not eliminate the local overdensity, allowing multiple cycles of oscillation. The 
result is a spectrum with multiple harmonics and thus multiple peaks. The 
existence and contribution of dark matter is only distinguishable from that 
of baryons alone with three or more peaks. Moreover, the magnitude of the 
third peak constrains the time of transition between a radiation and matter 
dominated universe. In particular, a finding that the second and third peaks 
were comparable in magnitude would suggest that dark matter dominated 
before recombination, which is a fundamental prediction of Big Bang cosmol- 
ogy. The magnitude of the third peak is also of interest for estimating the 
fraction of dark matter in the Universe. Astronomers have several methods 
for inferring the dark matter fraction (e.g., studying the rotation of galactic 
disks in the recent Universe), and it is vital to determine if these estimates 
are comparable to those produced by the physics of the early Universe. 

Finally, the pattern of CMB hot and cold spots we see on the sky corre- 
sponds to those photons just reaching us from the moment of recombination. 
(Recombination was relatively quick but not instantaneous, so there is some 
blurring of high spatial frequencies from the scatter of photons during that fi- 
nite period.) The contribution to this pattern from each acoustic mode maps 
to a spherical mode of fluctuations on the sky. The analysis then proceeds 
by decomposing the observed fluctuations into spherical modes and using 
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the contributions of these modes to understand the acoustic oscillations. We 
discuss this in the next subsection. 

1.2 The CMB Temperature Power Spectrum 

Our focus in this paper is inference about the CMB temperature power spec- 
trum and in particular the peaks in the spectrum. In this section, we describe 
the spectrum and some of the issues that arise in estimating it. Marinucci 
(2004, this issue) gives a more complete derivation upon which ours is based. 

Let T(0, i?) denote the temperature field as a function of colatitude (zero 
at the zenith) < 6 < n and longitude < $ < 2n. Let T denote the average 
of the temperature field over the sphere. 

Define the temperature fluctuation field by 



Note that Z is a random field with mean zero and is assumed to have finite 
second moment. We can expand Z in terms of a orthonormal basis on the 
sphere. The usual choice of basis is the set of spherical harmonics {Y^ m {9, ■&)}, 
for positive integers £ = 1,2,... and integers —£< m < £. (Here £ is called 
the multipole index, or loosely "multipole moment.") These are defined as 
follows: 



where the P^ m £ = 1,2,... and m = 0, . . . , £ are the associated Legendre 
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functions denned by 



Jm 



with Legendre polynomials 

2W dx i 



Pi(x) = ^^(x 2 ~iy. 



We can now write 

oo I 

Z(9,$)=J2 E <*l,mYt,m(e,#), (1) 



1=1 m= 

where, 



a e , m = [ [ Z{0,ti)Y tim {0,ti)smOffld&. (2) 
Jo Jo 

Since Z is a mean zero random field, the coefficients a^ m are random vari- 
ables. They have mean 0, variance 

Ci = E |a^ m | , 

and are uncorrected. The power spectrum is defined to be Cg as a function 
of t. 

Usually, it is assumed that Z is a Gaussian field (but see the paper by 
Marinucci in this issue) which implies that the a^ m have a Gaussian distribu- 
tion. If we were to observe Z without measurement error, we could estimate 

Ci by, say, 

1 1 

m=—£ 

and thus for large I we have Ce ~ Ce because we are averaging a large number 
of a| m . We call Cg, the realized spectrum. Another important implication of 
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equation (jSJ) is that even with perfect observations, we would not know the 
true power spectrum. Because our Universe is viewed as one realization of 
a stochastic process, Ct will in general differ from C#, especially for small £. 
This is known as the problem of cosmic variance. We return to this point in 
Section 

In practice, the data are subject to various sources of measurement error, 
blurring, and unobserved parts of the sky. For example,the Milky Way, which 
is relatively bright, obscures the deep sky along a wide band. The spherical 
harmonics are no longer orthogonal over what is left of the sphere, which 
induces correlation and bias into the estimated Ces. There are in addition a 
host of other complications in measuring Z. 

Our model, in vector form, is 

C = C + e, (4) 

where C is the observed spectrum and the noise vector e, with covariance ma- 
trix CKC T , incorporates the known sources of error, including measurement 
error. If there were no sky cut for the galaxy, K would be diagonal, but in 
practice, it incorporates the various known sources of error. In practice, the 
unknown C in the covariance matrix is replaced by a a pilot estimate, C . 
The choice of C° turns out to have surprisingly little effect on the results. We 
thus take the covariance matrix of e in equation (|3J) to be known and equal 
to E = C°K(C°) T 

Another issue is that the observations are actually derived from a convo- 
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lution of the C^s with £- dependent window functions; that is, the model is 
actually C = KC + e for some matrix K . But as Figure H3 shows, the rows 
of K are very nearly delta functions. (See Knox, 1999.) And in fact, incor- 
porating these window functions has negligible effect on our results, so we 
disregard them in what follows. 

2 Uniform Confidence Sets For Nonparamet- 
ric Regression 

Taking Yi = Ce and xi = £/L max , let f(xe) = Ce denote the true power 
spectrum at multipole index i. See Figure |U for the Ye from the WMAP data 
(Hinshaw et al. 2003). We can then rewrite equation © in the form of a 
nonparametric regression problem: 

Yi = f(xi) + eg, £ = L min , . . . , L max , (5) 

where e = (eL min , • • • , CL max ) is assumed Gaussian with known covariance ma- 
trix S as described earlier. This is only an approximation to the model ac- 
tually used, but we will not discuss the various practical complications here. 

Let aj denote the diagonal elements of S and n = L max — L min + 1 be the 
total number of observed multipoles. Henceforth, we will use % = I — L min + 1 
as an index. 

Our approach is to nonparametrically estimate the regression / and find 
a nonparametric 1 — a confidence ball B n for /. More precisely, we want B Tl 

liminf inf P(f E B n ) > 1 - a (6) 
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for some large function class T such as a Sobolev space. 

Once we have computed the confidence ball, we can construct a confidence 
interval for any functional T(f) of interest, such as the location of the first 
peak. If T is a set of such functionals and 

/„(T) = (mi„r(/), s r(/)) 

then we have that 

liminf inf P (T(f) G I n {T) for all T e T) > 1 - a. (7) 

Alternatively, we can construct the set of cosmological parameters that pro- 
duces spectra within the confidence ball, which gives a joint confidence set 
on these parameters. 

We use orthogonal series regression to estimate / and then construct 
a confidence ball via the Beran-Dumbgen pivot method (Beran 2000, and 
Beran and Diimbgen 1998), which was inspired by an idea in Stein (1981). 
Specifically, we expand / in the cosine basis / = YlJLo ^j'Pji where (f>o(x) = 1 
and (f)j(x) = \^2(f)(irjx) for j > 1. If / is fairly smooth, for example if / 
lies in a Sobolev space, then X^ >n /^ * s negligible and we can write f(x) 

1 " 

^r-VW (8) 

n 

i=i 

for < j < n. Note that vector Z is approximately Normal with mean fi 
and variance matrix UHU T / y/n, where U is the cosine basis transformation 
matrix. We define the monotone shrinkage estimator by 

% = \ Z i ( 9 ) 
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where 1 > Ai > • • ■ > \ n > are shrinkage coefficents. The estimate of / is 

n 

In this paper, we will use a special case of monotone shrinkage, called nested 
subset selection (NSS), in which Xj = 1 for j < J and Xj = for j > J. In 
this case, 

j 

f( x ) = ^2 Z $A X )- 
j=i 

The squared error loss as a function of A = (Ai, . . . , A n ) is 
Ln(X) = J (f(x) - f{x)fdx » Yfri - ^') 2 - 



The risk is 



R(X) = E f (f(x) - f{x)fdx » £ Xfl + f> - A,) 2 /i 

^ , = 1 n ,=1 



where a 2 = V(ej). The shrinkage parameter A is chosen to minimize the 
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Stein's unbiased risk estimate 



«w = E A ?f + E( 1 - A ,-) 2 (^-f) ■ 

7=1 7 = 1 V 7 + 



(10) 



3=1 3~- 

Beran and Diimbgen showed that R(X) is asymptotically, uniformally close 
to R(X) in either the monotone or NSS case. 

The Beran-Diimbgen method is based on the weak convergence of the 
"pivot process" B n (X) = y/n(L n (X) — R(X)) to a Normal (0, r 2 ) for some 
t 2 > 0. (The estimator for r 2 is given in the Appendix 3.) It follows that 
n , L n (X n ) - S n (X n ) 
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is an asymptotic 1 — a confidence set for the coefficients, where z a denotes 
the upper a quantile of a standard Normal and where jttj = /2j(A n ). Thus 

B n =^f(x) = J2fi J <P J (x): /iGD n | (11) 

is an asymptotic 1 — a confidence set for /. 

The approach to confidence sets that we use here is quite different than 

the more familiar confidence band approach in which one constructs bands of 

the form f(x) ± c\J\/ar(f(x)) for some c. The advantage of bands is that by 

plotting them, we get a simple visual impression of the uncertainty. However, 

there are some drawbacks to bands. In their most naive form, the constant 

c = z a /2, which does not account for the multiplicity over the xs. This can be 

fixed by using a larger constant, although the computation of the constant is, 

in some cases, nontrivial. See Sun and Loader (1994). Second, the available 

results about coverage appear to be pointwise rather than uniform over / e T 

although we suspect that the results can be strengthened to be asymptotically 

uniform. The third, and most serious problem, is that the function estimate 

/ is biased so the confidence interval is not centered properly, resulting in 

uncercoverage. Specifically, letting six) denote the standard error of / and 

m(x) = E f(x), we have that 

f(x) - f(x) = f(x) - m(x) m(x) - f(x) 
s(x) s(x) s(x) 

The first term typically satisfies a central limit theorem. The second term 
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does not tend to zero since optimal smoothing causes the bias m(x) — f(x) 
to be of the same order ). There have been some attempts to control 

this smoothing bias; see Ruppert, Wand and Carroll (2003) for a discussion. 

The confidence ball approach automatically deals with the smoothing 
bias, at least approximately. This is because the ball takes the object \ \f(x) — 
f(x)\\ 2 as its starting point, rather than f(x) —m{x) which is implicit in the 
band approach. The ball approach does have some bias, since / actually 
estimates f n (x) = Y^j=i rather than f(x) = resulting 

in a tail bias of y^°° , n u 2 -. However, this tail bias is small relative to the 

i— /J =71+1 r~J " 

smoothing bias. 

3 Dealing with Heteroskedastic Errors 

As Figure shows, the data for the CMB power spectrum are highly het- 
eroskedastic. The confidence set based on C? loss is a ball and thus gives 
equal weight to deviations in all direction. Because the the CMB variances 
are tiny for some is and huge for others, this symmetry is inappropriate. In 
parametric inference, confidence sets under heteroskedasticity are typically 
ellipses rather than balls, and we need to make a similar adjustment. We do 
this by constructing the confidence set under a loss function that gives more 
weight to points where the spectrum is measured precisely. In this section, 
we extend the Beran-Dumbgen method to such weighted loss functions. 
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We now replace the L 2 loss function with the following weighted loss: 



where we take w 2 (x) = l/a 2 (x). We expand both the unknown function and 
the weight function w 2 in the orthonormal basis. Hence, we write 



where <po, 4>x, ■ ■ ■ is the cosine basis on [0,1] defined above. 

The construction of B n requires a new central limit theorem and a mod- 
ified estimate of the asymptotic variance. We also replace the risk estimator 
in equation (fTU|) by the following, which can be shown to be unbiased for the 
new loss function: 



where D and D = I — D are diagonal matrices with l's in the first J and 
last n — J entries, B = XJYJJ T is the covariance of Z, and Wjk = w^Aj^ 
with wg being the £ th expansion coefficient of the function w 2 and 





R = Z T DWDZ + tiace(DWDB) - trzce(DWDB) 



(12) 
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The set B n is defined as in equation (jlljl but with the new estimate of risk. 
The estimated variance of the pivot, r 2 , is also different and is given in the 
appendix. 

4 Results 

We applied our method to the WMAP data to obtain a confidence set for 
the unknown spectrum /(£/L max ) = Ci. Figure |H1 compares the center of our 
confidence ball with the so-called "Concordance model" (Spergel et al. 2003). 
The Concordance model is the maximum likelihood estimator for a likelihood 
of the form 

^Conc(^) ^WMAP) ^LSS; Kyman) ^CBI ; ^Acbar) 
= £ W MAp(^) ^WMAp) ' 

^LSs(^) ^LSs) ' ^Lyman(^i ^Lyman) ' ^GBl(^) ^CBl) " -^-Acbar(^) ^Acbar) 1 (13) 

where the Ys are independent data sets from different experiments (WMAP: 
Bennett et al. 2003; LSS: Percival et al. 2003; Lyman: Croft et al. 2002, 
Gnedin and Hamilton 2002; CBI: Mason et al. 2003, Sievers et al. 2003, 
Pearson et al. 2003; Acbar: Kuo et al. 2001). In particular, Y WMAP is the 
data set we are using. The parametric fit from the WMAP data alone (see 
Figure E3 top right) is obtained by maximizing only the first component 

^WMAp(^) ^WMAp)' 

Note how well the nonparametric curve compares to the Concordance 
spectrum. The notable exceptions are in the very high-£ region around the 
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third peak and the \ow-£ region where the physical models curve upward 
sharply. We will argue that both the third peak and the rise in the spectrum at 
low is are by-products of the model and not the data. All of the cosmological 
models share both features. We are not suggesting that these features are 
incorrect, but we believe it is useful to separate effects driven by the data 
from those driven by the model. See Section |3J 

Once we construct the confidence ball, the next step is to use it to draw 
inferences. Because the ball is 900-dimensional in this case, it can seem daunt- 
ing to display results. Fortunately, our construction provides simultaneous 
coverage over all functionals of the unknown function, pre or post hoc. We 
thus explore the uncertainty by creating targeted probes of the ball. 

First, we explore the uncertainty in the location and height of the peaks. 
To do so, we searched through the confidence ball using local quadratic probe. 
Specifically, at each location £ , we defined a quadratic qu{£) with support on 
the interval [£q — A, £ + A] for fixed A = 51, centered at £q, and with height 
h. If / is the center of our confidence ellipse, we considered perturbations 
of the form / = f + qh- We varied h to find the largest and smallest values 
such that the resulting / is within the confidence ball and maintains three 
peaks over the I range [2, 900]. This results in confidence limits on the peaks 
as shown in Figure 

One striking result is the different precisions with which the first and sec- 
ond peaks are resolved. This is to be expected given the large variances near 
the second peak. In other words, the data alone give little information about 
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the second peak. (The third peak is even more uncertain.) The published 
results in the physics literature present the second peak with much lower 
uncertainty. We return to this issue in Section |5J 

Figure |H1 shows an example of a model-directed probe. Using the CMB- 
FAST software package (Seljak and Zaldarriaga 1996), we generated spectra 
in a one-dimensional family centered on the Concordance model parameters. 
The figure, which we call a ribbon plot, shows how the spectrum changes as 
the baryon fraction Q h is varied while keeping the total energy density fi Total 
fixed at 1. The light gray curves are in the ball; the black curves are not. The 
resulting interval for flbh 2 is [0.0169,0.0287]. To generate a valid confidence 
interval with such a probe we would need to search the entire 11-dimensional 
parameter space. 

The data are much noisier for high is, and we want to quantify how this 
propogates into local uncertainty about the function because this our ability 
to resolve the second and third peaks. A simple probe of the confidence 
set is useful for this purpose: finding how far a particular function in the 
confidence ball can be perturbed by localized deviations. For example, at 
each £, we examined the one- dimensional family of spectra fa = f + h ■ b, 
where b is a boxcar of fixed width and unit height centered at £. Figure El 
shows the maximum absolute height h that remains in the 95% ball relative 
to the height of the Concordance spectrum, for two different boxcar widths. 
At is where this curve is greater than 1, the data arguably contain little 
information about the height of the curve near that location. 
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The confidence ball is also useful for model checking. Figure El shows 
four different spectra along with the minimum value of 1 — a for which each 
spectrum is in the 1 — a confidence ball. The concordance spectrum is very 
close to the center. But the best fitting parametric model using only the 
WMAP data is at best in the 73% confidence ball. Cosmologists often use 
68% confidence levels, so this can be seen as a weak rejection of the best 
fitting model from the data. We also considered two extremal models that 
are in the 95% ball. These show that the data alone are consistent with 
eliminating the second and third peaks. While the cosmological models all 
predict these peaks - through the acoustic oscillations caused by dark matter 
- this suggests the benefits of more precise data, as from the Planck mission 
(Balbi et al. 2003). 

5 Discussion 
5.1 Findings 

Our most striking finding is that the center of our nonparametric confidence 
ball using the WMAP data alone lies very close to the Concordance model 
fit over the range where the data are not noise dominated. Recall from equa- 
tion (|13j) that the Concordance model incorporates four independent data 
sets, each with distinct likelihood forms. In contrast, the parametric model 
fit using only the WMAP data (with likelihood £wmap) lies barely in the 
73% confidence ball. Given that cosmologists often use 68% confidence in- 
tervals as their standard of evidence, this is tantamount to a rejection of the 
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cosmo logical model that underlies that parametric fit. 

This raises two points. First, it is remarkable that with a fully nonpara- 
metric method we have come very close to the Concordance model based 
on the WMAP data alone. Second, that we obtained basically the same 
spectrum as the Concordance model calls into question the accuracy of the 
WMAP-only likelihood £ WMAP . 

5.2 Methods 

We have presented a nonparametric method for analyzing the CMB spec- 
trum. Our techniques have wide applicability to regression problems beyond 
cosmology. By starting with a confidence ball, then probing the ball using 
functionals, one can address a variety of questions about the unknown func- 
tion while maintaining correct coverage, despite multiplicity and post-hoc 
selection. 

The method in this paper modifies the original Beran-Diimbgen construc- 
tion to account for heteroskedasticity. This modification yields a substantial 
reduction in the size of the confidence set. The resulting confidence set is also 
more useful in that it leads to tighter (looser) bounds in regions where the 
function is more (less) accurately measured. 

One advantage of our approach, is that it allows one to separate the 
information in the data from the information in a model. Although we did 
not pursue the full calculation here, we could intersect our confidence ball 
with the manifold of spectra from the parametric model as a way to combine 
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data and model. Specifically, we could use the cosmological model to generate 
spectra, but then test which spectra are consistent with the data by reference 
to our confidence ball. This does not rely on likelihood asymptotics which, 
as we discuss below, are suspect in this problem. Another advantage is that 
by extending this analysis to a constrained noparametric model (such as a 
three peak model) that contains the cosmological model, we can make the 
same inferences without being tied to the analytic form of the model. Our 
approach can then be used to check the model, make inferences under the 
model, and compare parametric to nonparametric inferences. 

We should point out that cosmologists obtain confidence intervals for pa- 
rameters in their (11-dimensional) model by integrating over the nusiance 
parameters and producing a marginal posterior. However, the likelihood is 
ill-behaved, under-identified, and degenerate. Moreover, in the physics liter- 
ature, there does not seem to be a clear appreciation of the fact that interval 
estimates obtained this way need not have correct frequentist coverage. 

There are several other advantages to our approach. If a parameter is 
under-identified this will show up automatically as a wide confidence interval. 
The intervals have correct asymptotic coverage and simultaneous validity 
over all parameters of interest. There is no need to integrate or profile the 
likelihood function. Finally, the asymptotic theory for (jHJ is insensitive to 
the fact that the standard asymptotics for the likelihood approach fail. 
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5.3 Inferential Foundations 

Interestingly, there seems to be some confusion about the validity of fre- 
quentist inference in cosmology. Since we have access to only one Universe 
- and thus cannot replicated it - some feel that it makes no sense to make 
frequentist inferences. This represents a common misunderstanding about 
frequentist inference in general and confidence intervals in particular. The 
frequency statements for confidence intervals refer to the procedure, not the 
target of the inference. Our method for constructing confidence balls traps 
the true function 95 percent of the time, even over a sequence of different, 
unrelated problems. There is no need to replicate the given experiment, or 
Universe. 

Complicating matters is the fact that the coefficients ag >m of the temper- 
ature field are random and unknown. To see the importance of this point, 
it is useful to make a finer distinction by defining the realized spectrum 
Cg = (l/(2£ + 1)) J2 e |a^ m | 2 , the "true spectrum" Cg — E (Cg) and the mea- 
sured spectrum Cg. Note that all our inferences have actually been directed 
at the realized spectrum. Some phyciststs find it disturbing to be making 
frequentist inferences about Cg since it is a realization of a random vari- 
able rather than a parameter in the usual sense. But this is no different than 
making inferences about a random effect in a standard random effects model. 

These confusions have led to an interesting movement towards Bayesian 
methods in cosmology. Of course, when used properly, Bayesian methods 
can be very effective. Currently however, the Bayesian interval estimates in 
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the physics literature seem questionable, being based on unfettered use of 
marginalizing over high-dimensional, degenerate likelihoods using flat priors 
chosen mainly for convenience. Indeed, an active area of research is finding 
corrections for such intervals to make them have correct coverage. Moreover, 
the potentially poor coverage of the Bayesian interval seems not to have been 
widely recognized in the Physics literature. 

Appendix 1. CMB Data. 

The CMB is composed of photons. The temperature of these photons (2.726 
Kelvin) means that the radiation will be at the microwave wavelengths. The 
light is collected via a dish (or reflector) and fed into either a (1) bolometer, 
which senses small temperature change as the photon hits the detector, or 
(2) a high performance transistor. In some cases (such as the aforementioned 
COBE experiment), the telescope is placed in orbit above the Earth. In other 
cases, the telescope is placed on a balloon and launched into the atmosphere. 
With careful attention paid to ground reflections, CMB telescopes can also 
be placed on the ground in regions where the atmosphere will contribute little 
contamination (like the South Pole). In all cases, there are a series of steps 
leading from the raw data collection to the final power spectrum estimation. 

The raw data are collected in a time stream. For each moment in time, 
the telescope records a temperature difference on the sky between two widely 
separated points. For example, one of these locations could be a fixed source 
of known temperature, thus allowing the temperature at the other point 
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to be calculated. However, the comparison location need not be fixed (or 
known) and the absolute temperatures can still be solved for iteratively, 
using previous measurements. Throughout this process, the pointing of the 
telescope needs to be accurately known (as a function of time), as well as the 
calibration of the temperatures, and also the instrument noise. 

Appendix 2. Cosmological Parameters 

The physics of the Universe on large scales is well described by a small 
set of cosmological parameters. We describe some of the most important 
parameters below. 

Normalized Hubble Constant h. The Hubble constant is the rate of the 
Universe's expansion. Specifically, H — ^ where a is the size of the Universe 
and d is the rate of change in a. "Constant" is a misnomer since this is a dy- 
namic quantity. The "Hubble constant" refers to the value of H as measured 
today (H ); this is often normalized and reported as h = i7 /100.0 with units 
km s" 1 Mpc -1 . 

Total Energy Density Jl To tai- ^Totai is the energy density of the Universe 
divided by the critical density of the Universe: p crit = 3c 2 Hq/8tiG at which 
the Universe would be geometrically flat. fi Total can be broken down into the 

SUm Of different Components, like ^baryons, ^dark matter, ^neutrinos- 

Cosmological Constant A. A is a constant that was added by Einstein into 
his general relativistic field equations to produce a static Universe. The con- 
stant was later dismissed as unnecessary after the discovery by Edwin Hubble 
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that the Universe is not static, but expanding. However, recent studies show 
strong evidence for a cosmological constant term. A acts as a negative pres- 
sure and thus might accelerate the expansion of the Universe. We often speak 
of the energy density component which is then included in the sum of 

^Total- 

Baryon Density Q b . This is the density component of baryonic matter 
in the Universe (e.g. protons, neutrons, etc). The fraction of matter density 
that is baryonic (over the total matter density of the Universe which includes 
baryons and non-baryonic dark matter) is often measured to be in the range: 
15%-20%. 

Dark Matter Density fl d . The majority of matter in the Universe is de- 
tected indirectly through it's gravitational effects. Since it cannot be seen 
or measured in the laboratory, it has been dubbed "dark matter" . Q d is the 
energy density component strictly due to dark matter. 

Neutrino Fraction, f u . The fraction of the neutrino density over the total 
matter density: f u = Q u /(Q b + fl d ). 

Optical Depth r. We know today that most of the hydrogen in the Uni- 
verse is ionized. So at some time after recombination, the Universe was re- 
ionized, t is the optical depth due to Thomson scattering up to a redshift 
of z < Zionization- ^ Zlon%zatlon) arn e dt where ot is the Thomson scattering 
cross-section and n e is the electron density. 

Spectral Index n s . The primordial matter density fluctuation spectrum is 
proportional to the scale size raised to the power n, the primordial spectral 
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index. On large enough scales, the the CMB temperature power spectrum's 
spectral index (n s ) is then close (or equal) to the primordial spectral index. 

The spectrum may be approximated numerically as a function of these 
parameters using the CMBfast software package (Seljak and Zaldarriaga 
1996). Figure |H1 shows spectra corresponding to a range of cosmological pa- 
rameter settings. For example, the location and amplitude of the first peak 
is related to the total energy density f2 To tai- The baryon fraction Q h and the 
spectral index n s drive the ratio of the amplitude of the first and second 
peaks. The ratio of the amplitudes of the second to third peaks depends on 
the density of matter (fl b + Q d + f u ), h, and n s . 

Appendix 3. Estimating r. 

Recall from Section |3] that the cosine basis is defined on [0,1] by 



4>o(x) = 1, 4>j{ x ) = V / 2cos(7T j 7'x) j > 0. 



Then, if j and k are distinct and positive, 



2 cos(7ijx) cos(7ikx) 



cos(7r(j + k)x) + cos(pi(j — k)x) 




If J > 



<jyj = 2cos 2 (7rjx) = cos(7r2jx) + 1 




'0- 
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Hence, 

fl if#O',M = 0} = 3 

lo if #{j,k,£ = 0} = 2 

I S jk 5 e + djrfok + 5 k i5 j if k, £ = 0} = 1 

We thus have that 

L(fJ) = J(f~f) 2 w 2 

j,k,e J 

= 52(Pj - h) (fa - A) J2 w * A w 

j,k I 
= ((3-p) T W{(3-P), 

where W jk = J2e w ^jkt- 

Let A = 1 — A and let D(x) denote the diagonal matrix with x along 
the diagonal. Write (3 = D(X)Z. Assume Z has a Normal^, B) distribution. 
Then, E/3 = D{\)(3 and since Cov (ft, fa) = XjX k B jk , Var 0) = D(X)BD(X). 
Then, 

EL = E0-(3) T W0-(3) 

= tmce(D(\)WD(\)B) + f3 T D(X)WD(X)(3. 

The latter quadratic form can be written as Ylj k PjPkXjX k Wj k . We obtain 
unbiased estimate L by replacing (5j(5 k with ZjZ k — Bj k . 

For convenience, let D denote D(X) and D denote D(X), The result is 

L = Z T DWDZ + trace(DWDB) - trace(DWDB). 
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It follows that, 

L-L = r Wf3-2Z T DWf3 + Z T DWDZ - 

Z T (I - D)W{I - D)Z - trace((W - DW - WD)B) 
= (3 T W(3 - 2Z T DW{3 + Z T (DW + WD - W)Z + trace((ZW + WD - W)B). 

Let A = DW + WD - W and C = 2DW{3. Then, 

Var(L-L) = Var {Z T AZ - Z T C) 

= Var (Z T AZ) + Var (Z T C) - 2Cov {Z T AZ, Z T C) 
= 2tmce(ABAB) + [3 T Q[3, 

where 

g/4 = ABA + WDBDW — 2ABDW 

= (DW + WD - W)B{WD + DW -W) + WDBDW - 2(DW + WD - W)BDW. 

Hence, plugging in unbiased estimates of the linear and quadratic forms in- 
volving (3, we get an estimate of the variance: 

? 2 = 2 tmce(ABAB) + Z T QZ - trace(Q5). (14) 

Acknowledgments 

The authors would like to thank the Executive Editor for his careful reading 
and helpful suggestions. 



31 



References 



Papers below with references to astro-ph numbers can be found on-line in 
multiple formats at http://xxx.lanl.gov/archive/astro-ph In addition, 
the following websites have useful introductory material about WMAP and 
the CMB: 

1. The WMAP home page 
http : //map . gsf c . nasa . gov/ , 

2. Wayne Hu's CMB Physics Pages 

http : / /background . uchicago . edu/ ~whu/physics/physics . html , 

3. Ned Wright's Cosmology Tutorial 

http : //www . astro . ucla . edu/~wr ight/cosmolog . htm , 

4. and an introductory review of CMB studies 

http : / /astron . berkeley . edu/~mwhite/ rosetta/ index . html 

Balbi, A., Baccigalupi, C, Perrotta, F., Matarrese, S., & Vit- 
torio, N. (2003). Probing Dark Energy with the Cosmic Microwave Back- 
ground: Projected Constraints from the Wilkinson Microwave Anisotropy 
Probe and Planck Astrophys. J., 588, L5-L8. 

Bennett et al. (2003). First Year Wilkinson Microwave Anisotropy Probe 
(WMAP) Observations: Preliminary Maps and Basic Results Astrophys. J. 

32 



Supp. 148, 1-27. 

Beran, R. (2000). REACT scatterplot smoothers: superefficiency through 
basis economy. J. Amer. Statist. Assoc. 95, 155-171. 

Beran, R. and Dumbgen, L. (1998). Modulation of estimators and con- 
fidence sets. Ann. Statist. 26, 1826-1856. 

Croft et al. (2002). Toward a Precise Measurement of Matter Clustering: 
Ly-alpha Forest Data at Redshifts 2-4, Astrophys. J., 581, 20-52 

Gnedin, N. Y. and Hamilton, A. J. S., (2002). Matter power spec- 
trum from the Lyman-alpha forest: myth or reality? Monthly Notices of the 
Royal Astronomical Society, 334, 107-116. 

Halverson, N.W and Leitch, E.M. and Pryke, C. and Kovac, 

J. AND CARLSTROM, J.E. AND HOLZAPFEL, W.L. AND DRAGOVAN, M. 
AND CARTWRIGHT, J.K. AND MASON, B.S. AND PADIN, S. AND PEAR- 
SON, T.J. and Shepherd, M.C. and Readhead, A.C.S. (2002). DASI 
First Results: A Measurement of the Cosmic Microwave Background Angular 
Power Spectrum. Astrophy. J. I, 568, 38-45. 

Hinshaw et AL. (2003). First- Year Wilkinson Microwave Anisotropy Probe 
(WMAP) Observations: The Angular Power Spectrum Astrophys. J. Supp., 
148, 135-159. 



33 



Hu, W. (1999). CMB Anisotropies: A Decadal Survey, in Birth and Evolution 
of the Universe: RESCEU 1999, |astro-ph/0002520| 

Hu, W. (2001). Physics of CMB Anisotropies, Warner Prize Lecture, Amer- 
ican Astronomical Society Meeting, San Diego, Jan. 2001, 
http : / /background . uchicago . edu/ ~whu/Presentations/warnerprint . pdf 



Hu, W. (2003). CMB Temperature and Polarization Anisotropy Fundamen- 
tals, Annals of Physics, 303, 203, |astro-ph/02 10 696 



Hu, W. and Dodelson, S. (2002). Cosmic Microwave Background Anisotropies, 
Ann. Rev. Astron. Astrophys., 40, 171, astro-ph/0110414 

Hu, W. AND Sugiyama, N. (1995). Toward Understanding CMB Anisotropies 
and Their Implications Phys. Rev. D, 51, 2199-2630. 

Knox, L. (1999). Cosmic Microwave Background Bandpower Windows Re- 
visited, Physic Rev. D, 60, issue 10, 103516, 5 pages. 

Kuo et al. (2003). High Resolution Observations of the CMB Power Spectrum 



with ACBAR, astro-ph/0212289 

Lee, A. T. et. al. (2002), A High Spatial Resolution Analysis of the 
MAXIMA- 1 Cosmic Microwave Background Anisotropy Data, Astrophys. J., 
561, L1-L6. 



34 



Marinucci, D. (2004). Testing For Non-Gaussianity On Cosmic Microwave 
Background Radiation: A Review Statistical Science, this issue. 

Mason et al. (2002). The Anisotropy of the Microwave Background to 
1 = 3500: Deep Field Observations with the Cosmic Background Imager, 
Astrophys. J., 591, 540-555. 

Netterfield, C. B. et. al. (2002). A measurement by BOOMERANG 
of multiple peaks in the angular power spectrum of the cosmic microwave 
background. Astrophys. J., 571, 604-614. 

Pearson et al. (2003). The Anisotropy of the Microwave Background to 1 = 
3500: Mosaic Observations with the Cosmic Background Imager, Astrophys. 
J., 591, 556-574. 

Percival et al. (2001). The 2dF Galaxy Redshift Survey: the power spec- 
trum and the matter content of the Universe, Monthly Notices of the Royal 
Astronomical Society, 327, 1297-1306. 

Ruppert, D., Wand, M. and Carroll, R. (2003). Semiparametric Regression. 
Cambridge University press. Cambridge. 

Schwarz, D.J. (2003). The first second of the Universe, Annalen Phys., 12, 
220-270, astro-ph/0303574| 

Seljak, U. and Zaldarriaga, M. (1996). A Line-of-Sight Integration 
Approach to Cosmic Microwave Background Anisotropies Astrophys. J., 469, 



35 



437-444. 

Sievers et AL. (2003). Cosmological Parameters from Cosmic Background 
Imager Observations and Comparisons with BOOMERANG, DASI, and MAX- 
IMA, Astrophys. J., 591, 599-622. 

Smoot, G. et al. (1992). Structure in the COBE differential microwave 
radiometer first-year maps. Astrophys. J., 396L, L1-L5. 

Spergel, D. et al. (2003). First- Year Wilkinson Microwave Anisotropy 
Probe (WMAP) Observations: Determination of Cosmological Parameters 
Astrophys. J. Supp., 148, 175-194. 

Stein, C. (1981). Estimation of the mean of a multivariate normal distri- 
bution. Ann. Statist. 9, 1135-1151. 

Sun, J. and Loader, C. R. (1994). Simultaneous Confidence Bands in 
Linear Regression and Smoothing, Annals of Statistics, 22, 1328-1345. 



36 




Figure 1. (top) The CMB as seen by the COBE satellite. The angular 
resolution of the satellite is about 10° and the various shades correspond 
to hot and cool spots with respect to the CMB blackbody temperature, 
(bottom) The CMB from the Wilkinson Microwave Anisotropy Probe. No- 
tice the high angular resolution. Also notice that the large-scale structures 
are apparent in both the COBE and the WMAP data. Image courtesy of 
the WMAP Science Team and available at the WMAP Mission website: 
http : / /map . gsf c . nasa . gov 
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Figure 2. Estimated CMB spectrum showing the three peaks of interest. 
The underlayed boxes give the data ranges and uncertainties from a variety 
of older CMB experiments, not including WMAP. From Hu (2001). 
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Figure 3. Bandpower windows from the matrix K centered on (left to right) 
1 = 100,200,300,400. 
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Figure 4. Y e as a function of t for the WMAP data. 
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Figure 5. Noise standard deviation as error bars on data (above) and as a 
function of I (below). 
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Figure 6. Center of our confidence ball (curve with sharp rise at right) and 
the power spectrum for the Concordance model (curve with three peaks). 
Note the striking agreement between the nonparametric fit and the paramet- 
ric fit. 
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Figure 7. Center of our 95% confidence ball with superimposed 95% intervals 
for the heights and widths of the first two peaks. 
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Figure 8. "Ribbon" probe of the confidence ball within the parametric model 
keeping f2 Total fixed at 1 and varying the Baryon fraction Vt h h 2 from 0.01225 
to 0.03675. 
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Figure 9. Height of local "box car probe" that is just in the 95% confidence 
ball, divided by the height of the Concordance spectrum, for two different 
box car widths. The horizontal line is at height 1. 
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Figure 10. CMB spectra: (top left) Concordance model fit, (top right) 
WMAP-only model fit, and (bottom) two extremal fits. 
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